% fields H for rho > a & rho < a

function [rho, phi, z] = func_H_inside(r, m, CONSTS)

    k0    = CONSTS.k0;
    a     = CONSTS.a;
    
    [coef, p1] = func_coef_of_fields(m, CONSTS);
    
    if (isnan(p1))
        rho = 0;
        phi = 0;
        z   = 0;
        return;
    end

    
    B1 = coef(1);
    B2 = coef(2);

    q1     = func_q1_from_p(p1, CONSTS);
    q2     = func_q2_from_p(p1, CONSTS);
    n1     = func_n1_from_p(p1, CONSTS);
    n2     = func_n2_from_p(p1, CONSTS);
    beta_1 = func_beta_1(p1, CONSTS);
    beta_2 = func_beta_2(p1, CONSTS);

    Q1 = k0*a*q1;
    Q2 = k0*a*q2;
    bessel_m_q1 = besselj(m, Q1*r);
    bessel_mp1_q1 = besselj(m+1, Q1*r);
    bessel_m_q2 = besselj(m, Q2*r);
    bessel_mp1_q2 = besselj(m+1, Q2*r);
    
    rho  = -1i*(B1*(p1*bessel_mp1_q1 - n1*beta_1*m*(bessel_m_q1./(Q1*r))) + ...
                         B2*(p1*bessel_mp1_q2 - n2*beta_2*m*(bessel_m_q2./(Q2*r))));
    
    phi  = -(B1*n1*(bessel_mp1_q1 - beta_1*m*(bessel_m_q1./(Q1*r))) + ...
                      B2*n2*(bessel_mp1_q2 - beta_2*m*(bessel_m_q2./(Q2*r))));
    
    z    = -(B1*q1*bessel_m_q1 + B2*q2*bessel_m_q2);
    
end
